%
% Data loading for paper




fid = fopen('Crati_EU33S3D.txt');
            CA = textscan(fid,'%s%f%f%f%f%f%f%f%f%f');
            fclose(fid);
            CAS = CA(1); 
            statName = CAS{:};
            ca8 = CA{8};    % test on 8-th column
            nrs = numel(ca8);
%             if sum(isnan(ca8)) == nrs
%                 CA(8:end) = []; % in this case, no error ellipses are provided
%             end
            ANC = CA(2:end);
            ncn = size(ANC,2);
            AN = nan(nrs,ncn);
            for k = 1:ncn 
                AN(:,k) = ANC{k};
            end
            A1 = [(1:nrs)' AN];
            nvar = ncn+1;

 
fileIm = 'CalGGGF.jpg';

[pathim,fim,~] = fileparts(fileIm);
fileVert = fullfile(pathim,[fim '.mat']);
IMS = struct('Image',fileIm,'Vertices',fileVert);


% ns = A1(:,1); 

e1  = A1(:,2); 
n1  = A1(:,3);
lv = length(e1);
ve1 = A1(:,5);
vn1 = A1(:,6);
vz  = A1(:,7);
vz(end) = NaN;  % KROT


%---- for the "arrow scale" visualization 
esca = e1;  
nsca = n1;
vesca = ve1; 
vnsca = vn1;

% A new point is generated in order to represent a velocity vector
% along the east axis (scale arrow)
std_esca = std(esca-mean(esca));
std_nsca = std(nsca-mean(nsca));
esca(lv+1)  = min(esca)-std_esca/2;
nsca(lv+1)  = min(nsca)-std_nsca/2;

LW = 2;

%--- data visualization ---
x1 = 566000;
x1 = GSS.DownLeftCorner(1);
x2 = 685200;
x2 = GSS.TopRightCorner(1);
y1 = 4315000;
y1 = GSS.DownLeftCorner(2);
y2 = 4405000;
y2 = GSS.TopRightCorner(2);

[xg,yg] = meshgrid(x1:200:x2,y1:200:y2);

vzg = griddata(e1,n1,vz,xg,yg,'linear');
vzgm = vzg;
% Ibunga = vzgm > 3;
% vzgm(Ibunga) = 5;
% Ibindi = (vzgm > 1.1)&(vzgm < 5);
% vzgm(Ibindi) = 1.1;
% bungabunga = [-2.25:0.25:1.25 2 3 4 5];

figure
hold on; 
% contourf(xg,yg,vzgm,bungabunga); shading interp; axis equal; colormap turbo;
surf(xg,yg,-99*ones(size(xg)),vzgm); view(2);
shading interp; axis equal; colormap turbo;
colorbar
ALA = axis;
plot(e1,n1,'xk','LineWidth',2,'HandleVisibility','off');

if lv < 1000
    for k = 1:lv
        
        se = sign(vesca(k)); sn = sign(vnsca(k));
        if isempty(statName)
            strStatk = num2str(ns(k));
        else
            strStatk = statName{k};
        end
        ht = text(e1(k)-se*std_nsca/15,n1(k)-sn*std_nsca/15,strStatk);
        if k < lv
            set(ht,'Color','k','FontSize',14);
        else
            set(ht,'Color','m','FontSize',14);  % KROT
        end
    end
end
title('VERTICAL VELOCITY (mm/yr)','FontSize',16);
set(gca,'FontSize',14);
xlabel('EAST (m)','FontSize',16); 
ylabel('NORTH (m)','FontSize',16);

gs_ReadBGImage(IMS,ALA);
hold off;
axis equal; 
axis([x1 x2 y1 y2]);
